%% EXERCISE 32: Cantilever beam - students's version

format short
addpath('ex32_toolbox_local')


%% Partition of unity option

%% LME options (Machine precision)
optLME.dim   = 2;
optLME.verb  = 0; %0:off
optLME.grad  = 1;           % Computation of the Gradient 0:OFF 1:ON
optLME.hess  = 0;           % Computation of the Hessian  0:OFF 1:ON
optLME.TolNR = 1.e-10;
optLME.knn   = 0;
optLME.Tol0  = 1.e-06;
optLME.gamma = 1.8;      % value of gamma to compute LME

% Nodal Integration options
%cubature order: 1 (1 GPts), 2 (3 GPts), 4 (6 GPts), 5 (7 GPts), 6 (12 GPts)
optGL.orderGL = 1;

%% Matrial Parameters
L  = 4;         % length
D  = 1;         % height
P  = -1000;     % external load
E  =  1e05;     % Young modulus
%nu = 0.499999;       % Poisson coefficient
parameters.L  = L;
parameters.D  = D;
parameters.E  = E;
parameters.nu = nu;
parameters.P  = P;
Lx = L;
Ly = 0.5 * D;

%% Convervence curves
%m      = 5;
m_id   = (2:m);      %index of the L2 error norm to perform the fit
h_n    = zeros(m,1);
L2_err = zeros(m,1);
%NChen  = [5 9 17 33 65];

for n=1:m
  fprintf(1,'n= %d  init',n);
  % Node Points
  Ny = ceil(NChen(n)/2);
  Nx = NChen(n)*4-3;

  h  = L /(Nx-1);
  optLME.spacing = h;

  % This line is important
  center  = 0.5*[Lx Ly];
  x_nodes = UniformGrid2D(Nx, Ny, Lx, Ly, center);

  nPts  = length(x_nodes);
  h_n(n) = 1/sqrt(2*NChen(n)*Nx);
  
  ids   = (1:nPts);
  id_bd = [ids(x_nodes(:,1)==0) ids(x_nodes(:,2)==0) ...
    ids(x_nodes(:,1)==Lx) ids(x_nodes(:,2)==Ly)];
  id_bd = unique(id_bd);

  % Sample points and gauss weights
  [x_samples w_samples conectivity gPts] = MakeGLSamples2D(x_nodes, optGL);
  
  t=cputime;
  fprintf(1,'\tgamma=%4.2f | order=%2d\n',optLME.gamma,optGL.orderGL);
  
  %% Nodes thermalization, nodes and samples adjacency structures are computed
  [beta_n range_n] = NodalThermalization(x_nodes, optLME);
  optLME.beta   = optLME.gamma/(h*h);
  optLME.beta_n = beta_n;
  optLME.range_n= range_n;
  
  %% The basis functions are computed
  % adjacency structure with the nearest neighbors nodes to each sample point
  % nodal shape parameter
  s_near = samplesAdjacency_fem (x_samples,conectivity,gPts);
  
  % Local-max entropy basis functions computation
%   optLME.s_near = s_near;
%   outLME = wrapper_lme(x_nodes,x_samples,optLME);
%   
%   p_samp  = outLME.p_samp;
%   dp_samp = outLME.dp_samp;
%   fprintf(1,'cputime LME   : %4.3f\n', cputime-t);
  
   % FEM basis function computation 
 outFEM = wrapper_fem (  x_nodes , x_samples , conectivity, gPts);
 p_samp  = outFEM.p_samp;
 dp_samp = outFEM.dp_samp;

  %% The displacement field is computed
  u_h = ex32_SolveSystem_sol(dp_samp,s_near,x_nodes,x_samples,w_samples,parameters,optLME);

  %% L2 norm
  u_num  = reshape(u_h,2,length(x_nodes))';
  u_num  = SamplesSolution(u_num,s_near,p_samp);
  planeState =1;
  u_sol  = ex32_AnalyticalSolution(x_samples,P,D,L,E,nu,planeState);
  L2_err(n)= sqrt(sum(sum((u_num-u_sol).^2,2).*w_samples))/sqrt(sum(sum(u_sol.^2,2).*w_samples))

end

%% Analysis of Results
% fit_num  = polyfit(log10(h_n(m_id)),log10(L2_err(m_id)),1);
% disp('- - - - - - - - - - - - - - -');
% fprintf ( 1, 'LME  : m=%7.4g  b=%7.4g\n', fit_num(1),  fit_num(2));
% 
% % Plots
% h_line=[0.01 0.1]*2;
% err_line=[0.0001 0.01]*2;
% fit_line=polyfit(log10(h_line),log10(err_line),1);
% fprintf ( 1, 'line : m=%7.4g  b=%7.4g\n', fit_line(1),fit_line(2));
% 
% figure(1);clf
% loglog(h_n,L2_err,'ro-',h_line, err_line,'k-', ...
%   'LineWidth',2,'MarkerSize',6)
% xlabel('h','fontsize',14,'fontweight','b')
% ylabel('Relative Error','fontsize',14,'fontweight','b')
% legend('LME','Location','Best','Orientation','horizontal')
% title(strcat('Relative Error in Norm L_2 -- \gamma_{LME}=',num2str(optLME.gamma)),'fontsize',16,'fontweight','b')
% set(gca,'XTick',0.02:0.02:0.1)
% %xlim([0.001 0.11])
% 
% 
% % Node and Sample Points
% figure(2)
% plot(x_nodes(:,1),x_nodes(:,2),'or',...
%   x_samples(:,1),x_samples(:,2),'.b');
% axis equal
% legend('Node points','Sample points');
% 
% 
% % Deformed and undeformed (numerical and analytical) configurations
% figure(3)
% u_sol = ex32_AnalyticalSolution(x_nodes,P,D,L,E,nu,planeState);
% u_sol = u_sol+x_nodes;
% u_num = reshape(u_h,2,length(x_nodes))';
% u_num = u_num+x_nodes;
% plot(x_nodes(:,1),x_nodes(:,2),'.r')
% axis equal
% hold on
% plot(u_num(:,1),u_num(:,2),'xb')
% hold on
% plot(u_sol(:,1),u_sol(:,2),'.g')
% hold on
% legend('Undeformed','Deformed (numerical)','Deformed (analytical)');
